function sh2(N)

  

  
  time_steps = [0:N-1]./(N-1)*100;

  function resB = B(z)
    if z == 0
      resB = 6;
    else
      resB = (1/3).*z.^(-1).*(1+z+z.^2).^(-2).*(2.*3.^(1/2).*pi+6.*z+3.^(1/2).*pi.*z+6.*z.^2+3.*z.^3+(-2).*3.^(1/2).*(2+z).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*(2+z).*log(1+z+z.^2));
    end
  end
  function resSk = Sk(z)
    if z == 0
      resSk = -4;
    else
      resSk = z.^3.*(1+z+z.^2).^(-1).*((-1)+4.*z.^(-3)+4.*z.^(-2)+2.*z.^(-1)+(1/3).*z.^(-4).*((-4)+z.*(1+z).*((-6)+z.*(2+z.*(2+z)))).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2)));
    end
  end
  function resM = M(z)
    resM = zeros(2,2);
    resM(1,1) = 2.*z;
    resM(2,2) = 1;
  end
  function res = deq(z, u)
    k   = u(1);

    res = zeros(2,1);
    kz = Sk(z) + 2.*(2+k);
    
    res(1) = kz;
    res(2) = 0;
  end
  uin = [0; 0];
  opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @M);
  [tout, uout] = ode15s(@deq, time_steps, uin, opts);
  plot(uout)
  uout


end
